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ABSTRACT 


This  paper  proposes  a stopping  rule  (Rule  2)  for  terminating  data 
collection  in  a queueing  simulation  experiment.  The  appeal  of  the  rule 
lies  in  the  fact  that  data  collected  in  this  way  can  be  used  to  compute 
interval  estimates  with  coverage  rates  that  compare  favorably  with 
theoretically  specified  rates.  The  rule  relies  on  a comparison  between 
a priori  information  on  the  activity  level  (traffic  intensity)  p and 
a corresponding  sample  estimator  computed  during  the  course  of  simulation. 
Experiments  with  simulations  of  the  M/M/c  queue  with  c * 1,2,4  and 
p = .7,. 8,. 9,. 95  were  conducted  to  evaluate  the  rule.  The  experiments 
used  a starting  rule  (Rule  1)  proposed  in  Adlakha  and  Fishman  (1979)  to 
reduce  bias  due  to  the  initial  conditions  and  also  used  the  autoregres- 
sive method  to  obtain  interval  estimates  of  the  steady-state  mean.  For 
p = .7,  .8,  .9,  the  coverage  rates  are  close  to  the  specified  theoreti- 
cal coverage  rates  and  are  higher  than  those  reported  in  the  literature 
for  other  methods  of  interval  estimation.  The  data  reveal  a degradation 
in  the  coverage  rate  for  increasing  values  of  activity  level.  For 
p = .95  the  coverage  rates  are  somewhat  lower  than  those  expected  theo- 
retically, indicating  room  for  some  improvement  in  technique.  The  sample 
sizes  used  to  obtain  the  coverage  rates  are  moderate  and  are  insensitive 
to  variation  in  the  number  of  servers  and  the  activity  level.  The  rule 
can  be  easily  generalized  to  a wider  class  of  queueing  simulations. 
Furthermore,  experiments  with  a fixed  truncation  starting  rule  and  a fixed 
sample  size  stopping  rule  clearly  demonstrate  the  superiority  of  using 
Rule  1 and  Rule  2 together.  This  is  very  encouraging,  for  It  indicates 


ii 


a procedure  now  exists  for  controlling  the  detrimental  effects  of  initial 
conditions  and  skewness  on  interval  estimation  in  queueing  simulations. 


1.  Introduction 
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This  paper  proposes  a stopping  rule  for  terminating  data  collection 
in  a queueing  simulation  experiment.  Data  collected  In  this  way  are  used 
to  compute  interval  estimates  (for  the  population  mean)  with  coverage 
rates  that  compare  favorably  with  theoretically  specified  rates.  This 
comparability  is  the  appeal  of  the  rule.  Although  most  analysts  of  output 
from  discrete  event  simulations  with  stochastic  Input  would  agree  that 
Interval  estimates  of  parameters  of  Interest  are  useful,  prior  attempts 
to  develop  procedures  for  computing  such  useful  estimates  have  met  with 
limited  success.  The  principal  impediment  has  been  the  Inability  of  pro- 
posed techniques  to  produce  empirical  coverage  rates  as  high  as  the  theo- 
retically specified  coverage  rates  In  sampling  experiments  designed  to 
evaluate  the  proposed  estimation  techniques.  A review  of  past  perform- 
ance appears  in  Adlakha  (1979). 

The  reasons  for  these  inadequacies  are  not  difficult  to  find.  Let 
6 be  a point  estimate  of  6 and  let  var(8)  be  a point  estimate  of 
var(e)  . Then  a 1-a  interval  estimate  of  fl  usually  takes  the  form 

0 ♦ oXar(e)  , 0 being  the  l-a/2  quantile  of  either  the  normal  or  a 
t distribution.  Now  the  writers,  as  well  as  others,  have  observed  two 
potential  sources  of  difficulty  with  this  Interval  estimate.  In  particular, 
one  customarily  treats  0 as  normally  distributed  when  It  often  exhibits 
positive  skewness  In  practice.  Also,  cov[0,  var(0)]  > 0 , Implying 


that  below  (above)  average  § occurs  with  below  (above)  average  var(©)  . 

These  observations  lead  to  the  Intuitive  conclusions:  Positive  skewness 
and  positive  correlation  Imply  that  0 falls  below  0 more  than  fifty 
percent  of  the  time  and  the  Interval  estimate  width  2Q/lrar(8)  Is  shorter 
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than  normal  theory  suggests  at  least  fifty  percent  of  the  time.  This 
asymmetry  of  the  center  of  the  Interval  estimate  and  Its  shorter-than- 
expected  width  lead  to  a lower  coverage  rate  In  practice  than  one  would 
expect  for  a normally  distributed  0 and  uncorrelated  § and  var(0)  . 

To  overcome  the  Inadequacy  of  past  suggestions,  this  paper  presents 
an  approach  to  interval  estimation  that  uses  ancillary  information  avail- 
able to  an  analyst  before,  during  and  after  execution  of  a simulation  of 
a queueing  system.  This  information  consists  of  arrival  and  service  rates 
or  their  corresponding  reciprocals,  mean  Interarrival  and  mean  service 
times.  The  principal  idea  is  to  compute  sample  quantities  such  as  sample 
mean  interarrival  and  service  times  and  compare  them  with  their  correspon- 
ding true  values.  When  the  deviations  between  sample  and  corresponding 
true  quantities  fall  within  specified  tolerances,  one  stops  or  terminates 
the  run  and  uses  the  data  collected  so  far  to  compute  an  interval  estimate. 
The  rationale  here  is  that  making  sample  quantities  representative  of  their 
corresponding  true  quantities  may  make  0 representative  of  6 . If  this 
behavior  materializes,  then  one  intuitively  expects  greater  success  in 
computing  adequate  interval  estimates  for  0 . One  additional  motivation 
for  this  use  of  <i  priori  information  is  that  it  is  available  for  a wide 
variety  of  queueing  simulations. 

The  paper  recommends  a e topping  rule  that  exploits  this  relationship 
between  the  theoretical  and  sample  activity  levels  in  a queueing  simulation. 
When  used  with  a storting  rule  described  in  Adlakha  and  Fishman  (1979)  for 
beginning  data  collection,  and  with  the  autoregressive  method  of  Interval 
estimation  (Fishman  1971),  the  proposed  stopping  rule  leads  to  coverage  rates 
comparable  to  theoretically  specified  ones  for  a wide  range  of  activity 
levels.  Also,  the  coverage  rates  are  higher  than  those  reported  in  the 
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literature  for  other  methods  of  Interval  estimation  and  require  smaller  sample 
sizes  (for  comparison  see  Law  1977,  Fishman  1978a).  We  demonstrate  this 
success  by  application  of  our  proposed  starting  rule  - stopping  rule  pro- 
cedure to  simulations  of  the  M/M/1,  M/M/2  and  M/M/4  queueing  simulations 
with  activity  levels  p - .7,  .8,  .9  and  .95  . Section  2 formulates  the 
problem  to  be  investigated,  describes  the  proposed  stopping  rule  and  how 
it  is  to  be  used  with  the  Starting  Rule  1 in  Adlakha  and  Fishman  (1979). 
Section  3 presents  results  on  coverage  rates  and  related  performance 
measures,  showing  the  favorable  performance  of  the  proposed  starting  rule  - 
stopping  rule.  It  also  includes  an  analysis  of  sampling  variation  in 
stopping  times  when  using  the  proposed  stopping  rule.  Section  4 compares 
the  proposed  starting  rule  - stopping  rule  procedure  with  several  poten- 
tially useful  alternatives  and  demonstrates  the  superior  performance  of 
our  proposal. 

2 . Problem  formulatjon 

Consider  a simulation  model  of  a queueing  system  with  c servers 
in  parallel,  independent  interarrival  times  with  mean  1/X  and  indepen- 
dent service  times  with  mean  1/m  . Let  denote  the  elapsed  time 
between  arrivals  of  jobs  1-1  and  1 and  let  S^  denote  the  service 
time  of  arrival  1 . Assume  that  the  simulation  begins  with  the  arrival 
of  job  1 to  an  empty  queue  and  c Idle  servers.  Let  denote  the 
»yet>-n  of  completion  1 where  system  time  denotes  waiting  time  plus 
service  time.  Assume  that  an  ultimate  objective  of  analysis  Is  to  compute 
an  interval  estimate  of  the  mean  system  time  from  a sample  record  of  system 
times.  Also  assume  that  qiven  the  choice  between  starting  data  collection 
In  an  undercongested  or  a congested  system,  one  prefers  the  congested  one. 

See  Adlakha  and  Fishman  (1979)  for  a discussion  of  this  point. 
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Starting  Data  Collection 

After  n completions  occur  during  a simulation  run,  one  can  esti- 

.1  n n 

mate  1/A  and  1/w  by  n J T4  and  n"  F S4  respectively. 

1-1  1 1-1  1 

These  estimates  are  unbiased  and  Independent  of  Initial  conditions.  Now 
the  activity  level  or  traffic  Intensity  for  a queueing  system  usually  Is 
defined  as 


p a arrival  rate/  number  of  servers  * service  rate 

* A/(c  • u>) 


for  which  one  estimate  Is 


Since  (5n  Is  usually  a biased  estimator  of  p an  alternative,  presum- 
ably more  desirable,  estimator  is 


A 


(1) 


since  E(pn)  is  in  principle  derivable  for  most  common  interrarival  and 
service  time  distributions.  For  example.  In  the  case  of  exponential 
Interarrival  and  service  times  E(f*n)  * pn/(n-l)  so  that 

Pn  ■ (n-l)Pn/*»  • 

In  the  present  paper,  data  collection  begins  with  system  time  T+l 


where 
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Starting  Rule  l:  T * min(n:  and  hold) 


and 


S,  = 


(lo|  n - p|  s 6)  , 

* £I.n-l>  • 
{ml  * n- 1 } , 


6 = specified  tolerance  > 0 , m = specified  integer  > 0 , I = {_(n-l )/mJ 

The  quantity  p.  is  the  local  sample  activity  level, 
i *n 


i 


i=ml+l 


Si 


n-ml-l 


I.n 


l 


n-ml 


i=ml+l 


Ti 


This  sample  activity  level  applies  for  the  M/M/c  queue  and  requires  an 
adjustment  in  its  correction  factor  (n-ml-l )/(n-mI ) for  distributions 
other  than  the  exponential.  Adlakha  and  Fishman  (1979)  describe  the 
detailed  benefits  of  this  rule. 

In  words.  Rule  1 requires  one  to  use  a sample  activity  level  based 
on  at  most  m past  completions.  The  quantity  I denotes  the  number  of 
times  one  needs  to  reset  pT  ; i.e.t  the  number  of  iterations  minus  one. 
A little  thought  shows  that 


pr(I  = i)  ■ (1  - qj'  q^, 


i = 0,1 .... 


where  qm  is  the  probability  of  success  on  a given  iteration.  Then  I 


has  a geometric  distribution  with  mean  (1  - q(n)/9m  and  variance 


(1  - q^)/q^  • Also,  the  mean  number  of  completions  E(T)  required  to 


• 1 
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meet  Starting  Rule  1 satisfies 


< E<T>  s m/Sn  ' 

Following  the  Adlakha  and  Fishman  (1979)  recommendations,  we  choose 
m * 5000  and  5 * .0001  . 


Selecting  a Stopping  Rule 

The  stopping  rule  that  we  develop  is  based  on  a comparison  of  the 
theoretical  activity  level  p and  the  sample  activity  level  pn  in 
(1)  computed  during  the  execution  of  a simulation  and  where  n now 
counts  the  number  of  completions  in  the  collected  data.  As  stated  earlier, 
we  anticipate  that  representing  p by  pn  will  induce  an  underlying  corre- 
lation structure  In  the  data  that  leads  to  a representative  sample  mean 
1 n 

Xn  = n £ XT+i  * where  xt+1 * XT+n  are  the  co11ected  system 
times  and  will  produce  an  Interval  estimate  or  confidence  interval  whose 

coverage  rate  agrees  with  the  corresponding  theoretical  coverage  rate. 

Although  the  notion  of  a stopping  rule  based  on  a comparison  of 
sample  and  corresponding  theoretical  activity  levels  is  easily  appreciated, 
several  Issues  play  an  Important  role  in  determining  the  explicit  form  a 
stopping  rule  should  assume.  These  issues  are: 

1.  What  should  the  tolerance  be  for  the  comparison  of 
Pn  and  p ? 

11.  Should  there  be  a directional  relationship  of  pn 

t0  Vi  ? 

111.  Should  there  be  a lower  bound  on  the  size  of  the 
sample  to  be  obtained  for  estimation? 

1v.  Should  there  be  an  upper  bound  on  the  sample  size  to 


avoid  excessive  data  collection? 
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We  address  each  issue  in  the  context  of  the  formulated  stopping  procedure: 

Begin  the  simulation  in  the  empty  and  idle  state. 

Apply  Starting  Rule  1 to  determine  the  first  system 
time  XT+1  to  collect. 

Collect  system  times  Xy+j , ....  Xy+N  where 

A 

N = min(n:  |pn  - p|  s 6*) 

and  pn  is  based  on  interarrival  and  service 
times  for  completions  T+l , ....  T+n  and  6* 
is  a specified  tolerance. 

Compute  a sample  mean  and  interval  estimate 
from  Xy^^  « ...  * . 

Tolerance  Criterion 

Step  3 implies  that  at  the  termination  of  data  collection 
PjY  C [p-6*,  p+6*]  . For  values  of  pn  C [p-6*,  p)  it  is  conceivable 
that  pn  < p-6*  for  all  n < N and  that  a congested  state  is  never 
experienced.  Alternatively,  when 

p < PN  ^ P + 6* 

it  is  possible  that  pn  > p + 6*  for  all  n < N and  that  the  system 
fails  to  experience  less  than  average  congestion.  Moreover,  the  ratio 
Pn/p  has  a positively  skewed  distribution  with  mean  1 and  is  indepen- 
dent of  p for  a fixed  6*  (see  Adlakha  1979).  This  indicates  that  as 


Stopping 
Rule  2 ^ 
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p Increases  pr(p_  < p - 6*)  Increases  faster  than  pr(p„  > p + 6*)  , 

Implying  that  It  Is  more  likely  that  p^  C To-6*,p)  than  p^  C [p,  p+6*]  . 

It  Is  reasonable  to  believe  that  If  the  choice  lies  between  over- 
estimating or  underestimating  the  congestion,  most  analysts  would  prefer 
to  overestimate  In  queueing  models.  This  Is  especially  Important  In  the 
case  of  systems  with  high  activity  level  p . To  favor  this  behavior  in 

1 

the  system,  we  modify  step  3 by  the  substitution:  j 

I 

N = minjn:  p s pn  < p t {*}  . j 

Directional  Criterion 

Once  pn  is  within  the  specified  interval  [p,  p+6*]  , it  indicates 
that  the  desired  congestion  level  has  been  reached.  At  this  stage  it 
would  be  preferable  to  continue  collecting  data  and  stop  only  when  there 
is  an  indication  that  the  system  is  leaving  the  congested  state,  i.e., 

A A I 

pn  < pn-l  ' We  exPect  this  to  happen  when  completion  n has  a shorter 
service  time  or  a longer  interarrival  time  than  the  average,  thereby 

I 

contributing  to  a below  average  system  time  for  the  nth  customer.  This 
will  also  produce  a decrease  in  pn  . Therefore,  we  include  an  additional 
condition  "pn  < Pn_i"  in  step  3 . We  anticipate  that  this  condition 
along  with  the  condition  "p  < pn  < p + 5*"  will  provide  a sample  record 
with  enough  observations  from  the  congested  system. 
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Loiter  Bound  on  the  Sample  Size 

A precise  estimate  of  var(Xn)  is  important  for  interval  es tuition. 
As  mentioned  earlier,  we  use  the  autoregressive  analysis  approach  to  esti- 
mate var(Xn)  . Based  on  experiments,  Fishman  (1971)  showed  that  this 
sample  variance  approximation  is  not  reasonable  when  the  sample  size 
n < 500  . To  improve  on  this  approximation,  we  restrict  stopping  time  to 
n > 1000  in  step  3 . 

Upper  Bound  on  Samf’le  Size 

Experience  in  testing  starting  rules  in  earlier  work  (Fishman  and 
Moore  1978)  revealed  that  in  the  absence  of  an  iterative  procedure  as  in 
Starting  Rule  1 excessively  large  starting  times  can  occur.  Step  3 
admits  a similar  possibility  for  stopping  times,  as  preliminary  experimen- 
tation confirmed.  To  favor  congested  systems  and  to  avoid  excessively 
large  stopping  times,  we  propose  an  iterative  stopping  rule  with  upper 
bound  of  m*  = 5000  , to  be  consistent  with  Starting  Rule  1 . The  upper 
bound  is  applied  as  follows:  If  the  stopping  rule  is  not  satisfied  until 
n = 5000  during  a simulation  run,  these  observations  are  discarded.  Another 
independent  replication  is  generated  from  the  empty  and  idle  state  and 
Starting  Rule  1 is  again  satisfied.  Then  the  complete  stopping  rule  is: 

N = min(n:  p < pn  s p+6*,  pn  < Pn_j.  n 2 1000,  m*  = 5000)  . 

A perusal  of  the  data  obtained  from  the  initial  experiments  indi- 
cates that  for  replications  where  pn_j  > p+6*,  the  empirical  coverage 
rates  are  comparable  to  the  theoretically  specified  coverage  rates.  A 
little  thought  shows  that  this  is  not  surprising.  It  is  reasonable  to 
believe  that  the  value  of  p^  might  play  a role  in  determining  the 


coverage  rate.  The  directional  criterion  p^  > ^ does  not  appear 
to  induce  sufficient  congestion  In  the  system  although  p s p^  s p+6*. 
because  It  is  possible  that  Pn  < P for  ns  N-2  . In  contrast,  the 
condition  p^ ^ > p+6*  favors  more  congestion  in  the  system. 


The  proposed  procedure  with  lower  and  upper  bounds  is  now: 

1.  Begin  the  simulation  in  the  empty  and  idle  state. 

2.  Apply  Starting  Rule  1 to  determine  the  first  system 
time  XT+^  to  collect. 

3.  n •*-  0 . 

4.  n «-  n + 1 . 


Stopping 
Rule  2 


5.  Collect  system  time  Xy+n  . 

^6.  If  ns  1000  go  to  4 . 

7.  Compute  pn  based  on  interarrival  and  system 
times  for  completions  T+l,  ....  T+n  . 

8.  If  p < p < p + 5*  and  P„  i > p + 6*  , N ■*-  n 

n n- 1 

and  compute  a sample  mean  and  interval  estimate  from 
! Xt+i  i . . . , and  stop . 

9.  If  n = m*  * 5000  , discard  Xj+^ , ...,  Xy+n  and  go 
v to  step  1 . 

10.  Go  to  step  4 . 


M 


Although  we  anticipate  that  the  upper  bound  m*  would  keep  the 
sample  size  from  becoming  excessive,  discarding  the  data  on  a replica- 
tion when  the  stopping  rule  Is  not  satisfied  may  affect  the  statistical 
reliability  of  the  estimate.  We  study  the  effect  of  this  upper  bound  In 
Section  3 . Preliminary  testing  Indicated  that  6*  ■ .01  would  be  a 
worthwhile  choice. 

We  study  Stopping  Rule  2 using  the  experimental  design: 

p ■ .7,  .8,  .9,  .95  , 
c * 1.  2,  4 . 

The  objective  here  Is  to  Investigate  the  performance  of  Stopping 
Rule  2 for  varying  number  of  servers  In  the  M/M/c  queue  over  a range 
of  values  of  p that  represent  moderate  to  high  levels  of  congestion. 

3.  Results  of  Experimentation 

To  study  Stopping  Rule  2 , we  performed  100  Independent  replica- 
tions on  each  queueing  model  c * 1,  2,  4 for  each  activity  level 
(p  = .7,  .8,  .9,  .95)  . On  each  replication,  the  collected  data  were: 

N*  « stopping  time.  Including  the  number  of  observations 
In  discarded  replications, 

XN  * sample  mean  system  time  based  on  N = N*  (mod  5000) 
observations, 

pn  = sample  activity  level  at  termination 

Rj  * sample  autocovariance  at  lag  j,  j * 0,1,..., 50  . 
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The  sample  mean  and  autocovariances  were  used  with  the  autoregressive 


method  to  produce  an  estimate  of  var(]fN)  and  90  and  95  percent 
Interval  estimates  for  mean  system  time. 

Coverage  Analysis 

Table  1 presents  the  empirical  coverage  rates  for  mean  system 
time.  The  data  show  that  the  coverage  rates  obtained  for  p s .9  are 
comparable  to  those  suggested  by  the  theory.  Notice  that  for  a given 
value  of  c , the  coverage  rate  tends  to  decrease  as  p increases. 

This  implies  that  higher  levels  of  utilization  In  the  system  degrade 
the  performance  of  Stopping  Rule  2 . For  a given  value  of  p , the 
performance  appears  to  be  unaffected  by  different  values  of  c . 

Since,  to  our  knowledge,  no  results  on  the  estimation  of  the  mean 
system  time  p for  the  M/M/c  queue  with  c > 1 and  p - .95  have  been 
reported  In  the  literature,  we  use  the  coverage  rate  for  the  M/M/1 
queue  with  p = .9  obtained  from  earlier  proposed  methods  as  our  basis 
of  comparison.  For  a sequential  method  developed  In  Fishman  (1971)  the 
90  percent  empirical  coverage  rates  obtained  were  between  66  percent 
and  79  percent.  In  an  empirical  testing  of  the  batch  means  method. 

Law  (1977)  obtained  a 90  percent  empirical  coverage  rate  of  86  for 
a fixed  sample  size  of  12800  . With  the  method  of  batch  means,  Fishman 
(1978a)  obtained  an  empirical  coverage  rate  of  90  for  the  95  percent 
confidence  Intervals  with  a fixed  sample  size  of  16384  . Therefore, 
we  can  conclude  that  Stopping  Rule  2 provides  a better  coverage  rate 
than  the  methods  cited  in  the  aforementioned  references.  The  significance 
of  this  achievement  is  magnified  by  the  fact  that  the  sample  sizes  used 
to  obtain  these  coverage  rates  are  much  smaller  than  those  used  by  other 
methods.  We  discuss  the  distribution  of  the  sample  size  shortly. 


In  general  one  would  like  sample  means  that  are  unbiased  estimators 
of  mean  system  time  u . However,  because  of  the  procedures  for  starting 
and  stopping  data  collection,  one  has  reason  to  believe  ii  priori  that 

XN  for  random  T and  N is  a biased  estimator  of  u . 

Let  Y,  denote  the  sample  mean  of  replication  1 for 


i * 1,  ...,  100  . Then 


I l 

? ' TO  V' 


gives  the  grand  sample  mean  over  all  replications  for  a given  (c.p) 
2-tuple  and  Y - p gives  an  estimate  of  the  bias  In  Y, , ....  Y1ftn  . 
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Then 

2 1 

s Y " T55^55 

is  an  estimate  of  var(Y)  . 

Table  2 shows  the  relative  bias  (Y  - u)/u  and  Table  3 shows 
the  standardized  deviate  (Y  - y)/»^y".  Although  the  significance  of 
the  bias  is  well  established  in  Table  3,  its  small  relative  impact  in 
Table  2 removes  any  substantive  concern  about  an  overriding  bias  com- 
ponent. 


100 

l (y.  - y y 

i-1 


Table  2 

Relative  Bias  of  the  Sample  Mean  System  Time4 
(Y  - u)/u 


p 

c N. 

.7 

00 

.9 

.95 

1 

.044 

.057 

.077 

2 

.048 

.076 

4 

.036 

.034 

.036 

.056 

4 The  mean  system  time  u can  be  computed  using  formulae 
in  Gross  and  Harris  (1974). 
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Table  3 

Standardized  Bias  of  the  Sample  Mean  System  Time 

(T  - u)/^ 


Xw  P 

c \ 

.7 

.8 

.9 

.95 

1 

4.409* 

2.792* 

1.874 

1.314 

2 

3.605* 

3.609* 

2.721* 

-1.287 

4 

4.879* 

3.243* 

1.411 

1.092 

* 

Significant  at  the  5%  level. 


Variance  Analysis 

When  computing  confidence  intervals  for  XN  , one  relies  upon  the 
accuracy  of  the  estimate  of  its  variance.  Here  we  examine  the  issue  of 

bias  in  the  estimation  of  the  variance  of  the  sample  mean  as  a possible 

"t* 

source  of  errorl  As  previously  discussed,  the  autoregressive  procedure 
was  employed  to  estimate  the  variance  of  the  sample  mean  in  a replica- 
tion.  Let  and  V.  denote  the  sample  mean  and  the  sample  variance 
of  the  sample  mean,  respectively,  on  the  1th  replication.  Furthermore, 
let 

A . 1 00  A 

V = ^i3l  Vi  * 


•j* 

Additional  comparisons  appear  in  Adlakha  (1979). 


, 100 

V * W £ <*i  - h‘ 


where 


, 100 

ree  I,.,  vi 


For  to  be  an  adequate  estimator  of  var(Y^)  , V should  be  close  to 

V since  V is  an  unbiased  estimator  of  var(Y^)  . 

A 

Table  4 presents  the  variance  estimates  V and  V . The  results 

A 

show  an  upward  bias  in  V for  p s .9  . This  pattern  is  somewhat  differ- 
ent for  p = .95  , where  in  two  out  of  three  cases  a downward  bias  in 

V is  suggested.  No  doubt,  this  underestimate  for  p - .95  is  responsible 
for  the  degradation  in  coverage  rate  in  Table  1 . Clearly,  work  remains  to 

a 

be  done  in  improving  V as  an  estimator  of  varfX^)  . 


Correlation 

Fishman  (1978a)  observed  that  a principal  source  of  error  in  the 
interval  estimation  of  correlated  data  is  the  high  correlation  between 
XN  and  var(XN)  . Table  5 provides  these  sample  correlations  for  the 
experiments  in  this  study.  The  data  confirm  that  and  var^) 
are  positively  correlated.  At  first  thought,  this  observation  Is  per- 
plexing, since  we  have  pointed  to  correlation's  being  the  culprit  that 
degrades  the  coverage  rate.  One  plausible  explanation  lies  in  the 


1/ 


llS 


Table  4 


Comparison  of  Estimates  of  var(7N) 


.7 

.8 

.9 

.95 

A 

V 

.107 

.580 

13.747 

64.700 

i V 

.060 

.393 

7.380 

124.557 

Ratio 

1.78 

1.48 

1.86 

0.51 

V 

.099 

.699 

11.292 

61.629 

2 V 

.053 

.345 

7.060 

59.568 

Ratio 

1.87 

2.03 

1.60 

1.03 

f 


significant  bias  In  X^  (see  Table  3).  This  small  but  perceptible 
upward  tendency  may  be  sufficient  to  make  the  correlation  In  Table  5 
work  in  our  favor.  That  Is  a tendency  to  overestimate  u with  a 
positive  correlation  between  and  var(7N)  leads  to  wider  Interval 
estimates  than  strict  normal  theory  would  suggest.  In  turn,  this  may 
contribute  to  the  improved  coverage  rates.  In  fact,  the  absence  of 
significant  bias  In  X^  for  p * .95  may  be  the  missing  Ingredient 

A 

needed  to  make  V more  representative  of  var^)  . 

Recall  that  Stopping  Rule  2 Imposes  an  upper  bound  m*  on  the 
sample  site  In  a replication.  To  study  the  affect  of  this  bound  on 
XN  , we  examine  the  sample  correlations  between  the  stopping  time  N* 
and  the  corresponding  sample  mean  y , N * N*(mod  5000)  in  Table  6 . 
The  correlations  show  no  significance  at  the  5 percent  level. 


Table  6 

Sample  Correlation3  Between  1L  and  N* 


p 

c 

.7 

.8 

.9 

.95 

1 

-.049 

.043 

.047 

2 

.106 

.192 

.068 

4 

-.049 

-.007 

.015 

.093 

a 

The  critical  value  with  100  observations  at  the 
5*  significance  level  Is  .195  . 
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Distribution  of  Stopping  Time 

Table  7 gives  the  sample  quantiles  of  the  stopping  time  N*  . 

The  lower  quantile  values  reflect  the  lower  bound  of  1000  in  Stopping 
Rule  2 . Although  the  quantiles  show  the  distribution  of  the  stopping 
time  to  be  positively  skewed,  the  95  percent  quantiles  range  from 
8722  to  17156  , indicating  small  probability  for  an  excessive  sample 
size. 

Notice  that  the  quantiles  appear  to  be  insensitive  to  p and  c . 

Since  the  distribution  of  pn/p  is  independent  of  c , the  insensiti- 
vity of  N*  to  c does  not  come  as  a surprise.  The  insensitivity  to 
p is  no  doubt  partially  due  to  the  facts  that  for  a fixed  6*  and  n , 

A 

as  p increases  the  pr(p  s pn  s p H*)  decreases,  whereas  pr(pp > p + 6*) 
increases.  The  effects  of  these  conditions  on  the  stopping  time  seem  to 
balance  each  other  for  an  increase  in  p . 

The  sample  mean,  standard  deviation  and  coefficient  of  variation  of 
the  stopping  time  N*  appear  in  Table  7 . The  mean  stopping  times  vary 
between  3586  and  5027  and  occur  between  the  60  and  70  percentiles. 

The  sample  coefficients  of  variation  0^*  are  generally  close  to  one. 

These  observations  suggest  that  although  the  distribution  of  N*  is  posi- 
tively skewed,  it  has  a short  tail  to  the  right.  Thus,  encountering  an 
excessive  stopping  time  is  a remote  possibility. 

4.  Comparison  with  Other  Rules 

This  section  compares  the  performance  of  the  proposed  starting-stopping 
rule  procedure  with  some  alternative  starting  and  stopping  rules  using  the 
experimental  design  which  we  describe  shortly.  We  use  the  coverage  rate  of 
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the  mean  system  time  as  the  criterion  for  comparisons. 

Let  Rj  and  R2  denote  Starting  Rule  1 and  Stopping  Rule  2 
respectively.  As  a variation  of  rule  R1  , consider  a fixed  truncation 
rule  that  puts  the  M/M/c  queue  into  the  steady  state  and  denote  this 
rule  as  F^  . Similarly,  as  an  alternative  for  rule  R2  , consider  a 
fixed  sample  size  rule  and  denote  it  as  . Rule  F$  uses  a sample 
size  equivalent  to  the  mean  stopping  time  of  rule  R2  for  given  c 
and  p . 

Four  combinations  of  starting  and  stopping  rules  are  considered: 

i.  (R-j , R2)  = Starting  Rule  1 with  Stopping  Rule  2 , 

ii.  (Fy,  R2)  = fixed  truncation  with  Stopping  Rule  2 , 

iii.  (R.| , F,. ) = Starting  Rule  1 with  fixed  sample  size, 

iv.  ( FT , F$)  = fixed  truncation  with  fixed  sample  size. 

Hereafter,  allusion  to  system  (A,B)  means  the  experiment  with  starting 
rule  A and  stopping  rule  B.  System  (•,  B)  refers  to  stopping  rule  B 
where  starting  rule  can  be  either  R^  or  Fy  and  system  (A,  •)  refers 
to  starting  rule  A with  stopping  rule  R2  or  F^  . 

Let  (Fq,  F<.)  denote  the  system  where  simulation  begins  in  the 
empty  and  idle  state  and  F$  is  the  stopping  rule  as  defined  earlier. 
Blomqvist  (1970)  showed  that,  provided  certain  conditions  are  satisfied, 
for  a large  sample  size  the  mean-square  error  of  sample  mean  as  an  esti- 
mator of  the  population  mean  is  minimized  if  one  starts  the  simulation 
from  the  empty  and  idle  state  and  no  observations  are  discarded.  However, 
for  a correlated  sample  record,  as  is  the  case  In  most  queueing  simula- 
tions, it  is  not  clear  that  minimizing  mean-square  error  would  necessarily 
provide  the  theoretically  specified  coverage  rate.  Therefore,  we  also 


\ 


i 


r — 
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compare  the  mean-square  errors  and  coverage  rates  obtained  by  using 
system  (Rp  R£)  with  those  obtained  from  system  (FQ,  F^)  . 

To  study  the  performance  of  the  aforementioned  systems,  we  consider 
the  experimental  design: 

p = .7,  .8,  .9,  .95  , 
c = 1,  2,  4, 

system  - (R^j,  R2)»  (Fy,  ^2^*  (fy»  • 

With  regard  to  truncation  rule  Fy  , an  exploratory  research 
showed  that  a truncation  of  1000  observations,  starting  from  the 
empty  and  idle  initial  state,  suffices  to  put  a system  in  the  steady 
state  for  the  M/M/c  queue  with  c = 1,  2,  4 and  p = .7,  .8,  .9  . 

This  implies  that  the  initial  conditions  hardly  affect  completion 
1001  for  these  queues.  However,  for  p = .95  , a truncation  of 
1000  observations  was  not  sufficient,  but  a truncation  of  2000 
observations  appeared  adequate. 

With  regard  to  stopping  rule  F<.  of  fixed  sample  size,  for  each 
given  value  of  p and  c we  take  a sample  of  size  N*  , the  corres- 
ponding mean  stopping  time  obtained  from  rule  . This  makes  rules 
R2  and  F$  comparable  in  the  sense  that  the  mean  number  of  observa- 
tions are  identical. 

Table  8 presents  the  results  for  100  replications  for  each  value 
of  p and  c with  systems  (Fy,  R2),  (R-|,  F^),  (Fy,  F^)  and  the  orig- 

inal results  for  (R.j,  R^)  . All  results  are  based  on  the  autoregressive 
analysis  in  Fishman  (1978b).  Three  tendencies  are  apparent.  Firstly,  the 
degradation  in  the  coverage  rate  for  all  the  systems  as  p increases. 


£ 


: 
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Table  8 

Coverage  Rates  for  Alternative  Systems 
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This  shows  that  high  utilization  in  the  M/M/c  queue  affects  the  performance 
uniformly  for  each  system.  Secondly,  the  superiority  in  the  performance 
of  system  (R^,  R2)  over  (R^ , Fs)  and  of  system  (F^,  R2)  over 
(FT,  Fs)  . This  indicates  that,  regardless  of  the  starting  criterion  used, 
stopping  rule  R2  yields  better  coverage  rates  than  the  fixed  sample  size 
F$  with  equivalent  expected  cost,  in  terms  of  number  of  observations. 
Thirdly,  there  is  no  apparent  distinction  in  the  performance  of  systems 


This  implies  that  the  performance  of  Starting  Rule  1 is  as  good  as  the 
performance  of  fixed  truncation  which  puts  the  system  in  the  steady  state. 

Two  subsequent  statistical  analyses  also  revealed  useful  insights. 

An  analysis  of  variance  together  with  multiple  comparison  procedures 
based  on  transformed  coverage  rates  showed  that  for  p s .9  , system 
(•,  R2)  was  significantly  better  than  system  (*,  F$)  . Moreover,  no 
significant  difference  existed  between  the  performance  of  systems  (R^  •) 
and  (Fy,  •)  . To  select  the  best  among  the  four  systems,  the  procedure 
in  Bechhofer  (1954)  was  used.  It  indicated  that  system  (R^,  R2)  was 
best  for  p = .7  and  .9  , with  at  least  probabilities  of  .562  and 
.956  respectively.  For  p = .8,  (F^.,  R^)  was  selected  as  best  with 
a correct  selection  probability  of  .458  . For  p * .95  , the  procedure 
was  indifferent  between  systems  (R^,  R2)  and  ( Fy . R2)  . The  prevalence 
with  which  the  several  statistical  analyses  focus  on  and  R2  enables 
us  to  conclude  that  (R^ , R^)  is  at  least  as  good  in  general  as  any  of  the 
other  systems.  However,  the  greatest  attraction  of  (Rp  R2)  comes  from 
the  fact  that  is  impractical  in  application  unless  one  knows  how  many 
observations  to  truncate  and  F^  is  impractical  unless  one  knows  how  many 
observations  to  collect  to  make  the  autoregressive  analysis  a useful  tool. 

As  mentioned  earlier,  the  Blomqvist  (1970)  result  encourages  one 
to  truncate  no  observations  to  minimize  mean-square  error.  To  test  the 
performance  of  the  zero  truncation  rule  Fq  , we  ran  the  design 
p = .7,  .8,  .9  and  .95  and  c = 1,  2 and  4 for  system  (Fq,  F^)  and 
computed  coverage  rates.  Table  10  compares  these  rates  with  those  for 

+See  Adlakha  (1979)  for  details. 
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(R.| , R2)  . The  superiority  of  (Rj , Rg)  is  clearly  established,  in 
spite  of  the  fact  that  computed  mean-square  errors  for  this  system  were 
larger  than  corresponding  mean-square  errors  for  (FQ,  F$)  . These 
observations  lead  to  the  conclusion  that  (FQ>  F$)  generates  a more 
accurate  point  estimate  in  terms  of  mean-square  error  but  (Rj,  R2) 
gives  a superior  interval  estimate  as  measured  by  coverage  rate.  In  our 
opinion  the  greater  precision  of  (Fq,  F^)  for  the  point  estimate  is  of 
little  value  if  we  have  no  adequate  way  to  estimate  the  precision.  We 
prefer  to  pay  the  penalty  of  a less  accurate  point  estimate  whose  preci- 
sion we  can  successfully  evaluate. 


Table  9 

Coverage  Rates  for  Systems  (FQ,  Fj)  and  (Rj , R?) 


90%  Coverage 


95%  Coverage 
P 
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5.  Conclusions 

The  objective  In  this  study  has  been  to  devise  generally  applicable 
methods  to  Improve  Interval  estimation  In  a queueing  simulation.  On  the 
basis  of  accumulated  empirical  evidence  presented  In  this  paper,  we  continue 
to  recommend 

Rule  1 : 

T * minjn:  |pj  n - p|  s 6,  Pj  n > pj  n_,,  ml  * n-1,  6 - .0001, 

m ■ 5000} 

as  a starting  rule  (see  Adlakha  and  Fishman  1979).  Also,  we  recommend 
Rule  2: 

N = minjn:  p s p>n  < p+6*.  pnJ  > p+«*,  n 2 1000,  m*  * 5000,  * .01} 

as  a stopping  rule. 

The  rules  give  satisfactory  performance  for  multlserver  queues 
with  moderate-to-high  levels  of  congestion.  Both  rules  use  a moderate 

number  of  observations,  are  easy  to  understand  and  to  Implement. 

Furthermore,  experiments  with  a fixed  truncation  starting  rule  and  a 

fixed  sample  size  stopping  rule  clearly  demonstrate  the  superiority 
of  using  Rule  1 and  Rule  2 together.  This  Is  very  encouraging,  for  it 
indicates  a procedure  now  exists  for  controlling  the  detrimental  effects 
of  initial  conditions  and  skewness  on  interval  estimation  In  queueing 
simulations. 
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tlons  of  the  H/M/c  queue  with  c ■ 1,2,4  and  p ■ .7,. 8,. 9,. 95  were  conducted 
to  evaluate  the  rule.  The  experiments  used  a starting  rule  (Rule  1)  to  reduce 
bias  due  to  the  Initial  conditions  and  also  used  the  autoregressive  method  to 
obtain  Interval  estimates  of  the  steady-state  mean.  For  p « .7,. 8, .9  , the 
coverage  rates  are  close  to  the  specified  theoretical  coverage  rates  and  are 
higher  than  those  reported  In  the  literature  for  other  methods  of  Interval 
estimation.  The  data  reveal  a degradation  In  the  coverage  rate  for  Increasing 
values  of  activity  level.  For  p ■ .95  the  coverage  rates  are  somewhat  lower 
than  those  expected  theoretically.  Indicating  room  for  some  Improvement  In 
technique.  The  sample  sizes  used  to  obtain  the  coverage  rates  are  moderate 
and  are  Insensitive  to  variation  In  the  number  of  servers  and  the  activity  leve 
The  rule  can  easily  be  generalized  to  a wider  class  of  queueing  simulations. 
Furthermore,  experiments  with  a fixed  truncation  starting  rule  and  a fixed 
sample  size  stopping  rule  clearly  demonstrate  the  superiority  of  using  Rule  1 
and  Rule  2 together.  This  Is  very  encouraging,  for  It  Indicates  a procedure 
now  exists  for  controlling  the  detrimental  effects  of  Initial  conditions  and 
skewness  on  Interval  estimation  In  queueing  simulations. 


UNCLASSIFIED 

SICUNITV  CLASSIFICATION  OF  THIS  FAStfOM*  Dm*  •»<•*•<> 


